home *** CD-ROM | disk | FTP | other *** search
/ The World of Computer Software / The World of Computer Software.iso / xlmath10.zip / XLMATH.C < prev    next >
C/C++ Source or Header  |  1992-03-14  |  19KB  |  663 lines
  1. /*********************************************************************
  2.  XLMATH.C
  3.  ========
  4.  This module demonstrates the basic structure of a Dynamic Link Library
  5.  that can be used with Microsoft Excel.
  6.  
  7.  Author:    Roy Kari
  8.              Dept. of Chemistry
  9.             Laurentian University
  10.             Sudbury, Ont.
  11.             Canada        P3E 2C6
  12.             (705) 675-1151
  13.             Internet: "ROY@NICKEL.LAURENTIAN.CA"
  14.  
  15.  ********************************************************************/
  16.  
  17. /* --------------------------< Include files >--------------------- */
  18.  
  19. #include <windows.h>
  20. #include <float.h>
  21. #include <math.h>
  22. #include <stdlib.h>
  23. #include "xlmopti.h"
  24. #include "xlmath.h"
  25. #include "xlmutil.h"
  26. #include "xlmcurve.h"
  27.  
  28. /* --------------------<forward references>------------------------ */
  29.  
  30. static WORD PASCAL BinBucket(LPFP lpValues, double dMin, double dMax);
  31. static int PASCAL Jacobi(LPLPREAL hmat, LPLPREAL eivec,
  32.                             LPREAL eival, int n);
  33. static void PASCAL SortVectors(LPLPREAL eivec, LPREAL eival, int n);
  34.  
  35. /* ----------------<Frequency utility routine>---------------------- */
  36.  
  37. /*********************************************************************
  38.     BinBucket()
  39.     ==============
  40.     This is an EXCEL math routine that calculates the number of occurences
  41.     greater than the lower bound and less than or equal to an upperbound.
  42.     This routine is used internally by Frequency() but can also be called
  43.     from EXCEL.
  44.  ********************************************************************/
  45. static WORD PASCAL BinBucket(LPFP lpValues, double dMin, double dMax)
  46. {
  47.  
  48.     WORD wValueRow;
  49.     WORD wBinCount = 0;
  50.  
  51.     // no memory allocated so don't bother checking initialization
  52.     for (wValueRow = 0; wValueRow < lpValues->wRows; ++wValueRow)
  53.     {
  54.         if (lpValues->Data[wValueRow] > dMin &&
  55.             lpValues->Data[wValueRow] <= dMax)
  56.         {
  57.             ++wBinCount;
  58.         }
  59.     }
  60.     return wBinCount;
  61. }
  62.  
  63. /*********************************************************************
  64.     Frequency()
  65.     ==============
  66.  
  67.     An Excel DLL math routine to calculate the freqency of occurences.
  68.     Each value in the bin block represents all values from it down to
  69.     the previous value. The first bin value represents the number of values
  70.     less than or equal to itself. The last bin value represents the number
  71.     of values greater than itself. The number of bin values returned is
  72.     always one greater than the number of bin values.
  73.  ********************************************************************/
  74. LPFP PASCAL FAR _export Frequency(LPFP lpValues, LPFP lpBrackets)
  75. {
  76.     WORD wValueRow, wBracketRow;
  77.     LPFP lpXLKInt;
  78.  
  79.     // Init XL transfer array 1 row greater
  80.     if ( (lpXLKInt =InitRetBuff(lpBrackets->wRows+1, lpBrackets->wCols)) == NULL)
  81.         return (LPFP) NULL;
  82.  
  83.     // zero the bin counts
  84.     for (wBracketRow = 0; wBracketRow <= lpBrackets->wRows; ++wBracketRow)
  85.     {
  86.         lpXLKInt->Data[wBracketRow] = 0;
  87.     }
  88.  
  89.     // check first and last limits
  90.     for (wValueRow = 0; wValueRow < lpValues->wRows; ++wValueRow)
  91.     {
  92.         // value is <= first bin?
  93.         if ( lpValues->Data[wValueRow] <= lpBrackets->Data[0] )
  94.         {
  95.             ++lpXLKInt->Data[0];
  96.         }
  97.         // value is > last bin?
  98.         else
  99.             if ( lpValues->Data[wValueRow] >
  100.                     lpBrackets->Data[lpBrackets->wRows-1])
  101.             {
  102.                 ++lpXLKInt->Data[lpBrackets->wRows];
  103.             }
  104.     }
  105.  
  106.     // check remaining bins if value > bin[-1] and <= bin[0]
  107.     for (wBracketRow = 1; wBracketRow < lpBrackets->wRows; ++wBracketRow)
  108.     {
  109.         lpXLKInt->Data[wBracketRow] = (double)BinBucket(lpValues,
  110.             lpBrackets->Data[wBracketRow-1],lpBrackets->Data[wBracketRow]);
  111.     }
  112.     return  ((LPFP)lpXLKInt);    // return as a long pointer
  113. }
  114.  
  115. /* ----------------<Diagonalize utility routines>------------------ */
  116.  
  117. /********************************************************************
  118.  Jacobi()
  119.  ========
  120.  This routine performs a diagonalization of a real symmetric matrix.
  121.  Only the the top half need be defined.  The routine is adapted from
  122.  a FORTRAN routine of the same name.  It is long and probably
  123.  should be broken into manageable portions but since it is
  124.  the biggest single time user in this application, I have opted
  125.  for speed over clarity.
  126.  *******************************************************************/
  127.  
  128. /* -------------------------< Defines >---------------------------- */
  129.  
  130. #define SQRTWO  0.707106781
  131. #define ACCRCY  DBL_EPSILON*1.e02   /* accuracy of diag        */
  132. #define MACHACC DBL_EPSILON            /* machine precision    */
  133.  
  134. static int PASCAL Jacobi(LPLPREAL hmat, LPLPREAL eivec,
  135.         LPREAL eival, int n)
  136.     {
  137.     int i, j;                              /* loop counters          */
  138.     int converged;                  /* convergence control    */
  139.     int jcol, irow;                 /* loop counters          */
  140.     int it_cnt;                     /* iteration counter      */
  141.  
  142.     double htop, dstop, thrsh, d;   /* stopping criteria etc   */
  143.     double fd, avgf;                /* temporary variable      */
  144.     double hii, hjj, hij;           /* elements of hmat        */
  145.     double c, s, t, u, dsu, dcu;    /* cos, sin, tan etc       */
  146.  
  147.     it_cnt = 0;
  148.  
  149.     /* set eigenvector matrix to unit matrix */
  150.     for( i = 0; i < n; i++)
  151.         {
  152.     /* set eivec to unit matrix      */
  153.     for (j = 0; j < n; j++)
  154.         {
  155.         eivec[i][j] = ( i==j ) ? 1.0: 0.0;
  156.         }
  157.         }
  158.  
  159.     /* find max element of abs(hmat) */
  160.     htop = 0.0;
  161.     for(i = 0; i < n; i++)
  162.         {
  163.     for (j = i; j < n; j++)
  164.         {
  165.         htop = (htop > fabs(hmat[i][j])) ? htop : fabs(hmat[i][j]);
  166.         }
  167.         }
  168.  
  169.     /* if only one element then done */
  170.     if (n == 1)
  171.         {
  172.     eival[0] = hmat[0][0];
  173.     return(0);
  174.         }
  175.  
  176.      /*  find the norm of the function and calculate stopping criteria */
  177.     avgf = (n-1)*n*0.55;
  178.     d = 0.0;
  179.     for (i = 1; i < n; i++)
  180.         {
  181.     for(j = i;j < n; j++)
  182.         {
  183.         s=hmat[i-1][j]/htop;
  184.         d = s*s + d;
  185.         }
  186.         }
  187.     dstop = d*ACCRCY;
  188.     thrsh = sqrt(d/avgf)*htop;
  189.  
  190.     while (TRUE)
  191.         {
  192.     /* MAIN LOOP */
  193.     converged = TRUE;
  194.  
  195.     /* loop over columns */
  196.     for (jcol=1; jcol<n; jcol++)
  197.         {
  198.  
  199.         /* loop over rows   */
  200.         for (irow=0; irow<jcol; irow++)
  201.             {
  202.  
  203.         hij = hmat[irow][jcol];
  204.         if(fabs(hij) > thrsh )
  205.             {
  206.             /* executed for off-diagonal element > threshold */
  207.             hii = hmat[irow][irow];
  208.             hjj = hmat[jcol][jcol];
  209.             s = hjj - hii;
  210.  
  211.             if (fabs(hij) > MACHACC*fabs(s))
  212.                 {
  213.             /*  executed for rotations greater then rounding error
  214.                 and any element nulls convergence
  215.             */
  216.             converged = FALSE;
  217.             if((double)(MACHACC*0.1)*fabs(hij) >= fabs(s))
  218.                 {
  219.                 /*
  220.                              *  If the rotation is close to 45 degrees,
  221.                              *  sin and cos to 1/root2.
  222.                              *  else calculate both sin & cos
  223.                              */
  224.                 s = SQRTWO;
  225.                 c = s;
  226.                 }
  227.             else
  228.                 {
  229.                 t = hij/s;
  230.                 s = 0.25/sqrt(0.25 + t*t);
  231.                 c = sqrt( s+0.5);
  232.                 s = 2.*t*s/c;
  233.                 }
  234.                /*
  235.                          *  calculate new elements of hmat
  236.                          */
  237.             for (i=0; i<=irow; i++)
  238.                 {
  239.                 t = hmat[i][irow];
  240.                 u = hmat[i][jcol];
  241.                 dsu = s*u;
  242.                 hmat[i][irow] = c*t-dsu;
  243.                 dcu = c*u;
  244.                 hmat[i][jcol] = s*t + dcu;
  245.                 }
  246.  
  247.             if( (irow+2) <= jcol )
  248.                 {
  249.                 for (i = irow+2; i <= jcol; i++)
  250.                     {
  251.                 t = hmat[i-1][jcol];
  252.                 u = hmat[irow][i-1];
  253.                 dsu = c*t;
  254.                 hmat[i-1][jcol] = s*u + dsu;
  255.                 dcu = s*t;
  256.                 hmat[irow][i-1] = c*u - dcu;
  257.                     }
  258.                 }
  259.             dsu = c*hjj;
  260.             hmat[jcol][jcol] = s*hij + dsu;
  261.             dcu = s*hjj;
  262.             fd = c*hij - dcu;
  263.             dsu = s*fd;
  264.             hmat[irow][irow] = c* hmat[irow][irow] - dsu;
  265.  
  266.             for (j=jcol; j<n; j++)
  267.                 {
  268.                 t = hmat[irow][j];
  269.                 u = hmat[jcol][j];
  270.                 dsu = s*u;
  271.                 dcu = c*u;
  272.                 hmat[irow][j] = c*t - dsu;
  273.                 hmat[jcol][j] = s*t + dcu;
  274.                 }
  275.             /*
  276.                          * begin calculation of eigenvectors
  277.                         */
  278.             for (i=0; i<n; i++)
  279.                 {
  280.                 t = eivec[i][irow];
  281.                 dsu = eivec[i][jcol]*s;
  282.                 eivec[i][irow] = c*t - dsu;
  283.                 dcu = eivec[i][jcol] * c;
  284.                 eivec[i][jcol] = s*t + dcu;
  285.                 }
  286.                      /* recalculate the norm d and compare with dstop */
  287.             it_cnt += 1;
  288.             s = hij/htop;
  289.             d = d-s*s;
  290.             if( d < dstop )
  291.                 {
  292.                 d = 0.0;
  293.                 for (i = 1;i < n; i++)
  294.                     {
  295.                 for(j = i;j < n; j++)
  296.                     {
  297.                     s=hmat[i-1][j]/htop;
  298.                     d = s*s + d;
  299.                     }
  300.                     }
  301.                 dstop = d*ACCRCY;
  302.                 }
  303.             thrsh = sqrt(d/avgf)*htop;
  304.                 }           /*  end rounding if    */
  305.             }               /*  end threshold if   */
  306.  
  307.             }                   /*  end row loop       */
  308.         }                       /*  end column loop    */
  309.     if (converged)
  310.         break;              /*  end infinite while */
  311.         }
  312.  
  313.     /* copy diagonal elements to eival */
  314.     for (i = 0; i < n; i++)
  315.         {
  316.     eival[i] = hmat[i][i];
  317.         }
  318.  
  319.     /* get in beta order */
  320.     SortVectors(eivec, eival, n);
  321.     return(it_cnt);
  322.     }
  323. /********************************************************************
  324.  *
  325.  * Name         SortVectors -- sort  a matrix
  326.  *
  327.  * Synopsis     SortVectors(eivec, vec, n);
  328.  *
  329.  *              double eivec[][] - eigenvector matrix
  330.  *                                    (sorted to correpond to eivec)
  331.  *              double vec[]     - eigenvalue matrix (determines sort order)
  332.  *              int n            - no. of atoms
  333.  *
  334.  * Description  (Shell) sort matrix eivec[][]  in increasing order of
  335.  *              values in vec[].
  336.  *
  337.  * Returns      eivec[][] and eival[]
  338.  *
  339.  ********************************************************************/
  340. static void PASCAL SortVectors(LPLPREAL eivec, LPREAL eival, int n)
  341. {
  342.     int gap;
  343.     register int i, j, k;
  344.     double temp;
  345.  
  346.     for (gap = n/2; gap > 0; gap /=2) {
  347.     for ( i=gap; i<n; i++) {
  348.         for ( j=i-gap; j>=0 && eival[j] < eival[j+gap]; j-=gap)  {
  349.  
  350.         temp = eival[j];            /* exchange vec values */
  351.         eival[j] = eival[j+gap];
  352.         eival[j+gap] = temp;
  353.  
  354.         for (k=0; k<n; k++) {   /* exchange corresponding eivec */
  355.             temp = eivec[k][j]; /* columns                    */
  356.             eivec[k][j] = eivec[k][j+gap];
  357.             eivec[k][j+gap] = temp;
  358.         }
  359.         }
  360.     }
  361.     }
  362. }
  363.  
  364. /********************************************************************
  365.  Diagonalize()
  366.  =============
  367.  Diagonalize a real symmetric matrix. Eigenvalues are returned
  368.  *******************************************************************/
  369. LPFP PASCAL FAR _export Diagonalize(LPFP lpHmat)
  370. {
  371.     LPFP lpXLKInt;
  372.     LPLPREAL lplpHmat;        // pointers to EXCEL allocated data array
  373.     LPLPREAL lplpEivec;        // handle to pointers
  374.     LPREAL lpEival;            // pointer to eigenvalue array
  375.  
  376.     WORD id;                // index
  377.     int nIter;                // no iterations in Jacobi
  378.     WORD wRows = lpHmat->wRows;    // saves typing
  379.     BOOL Success = TRUE;
  380.  
  381.     // error if not square array
  382.     if (wRows != lpHmat->wCols)
  383.     {
  384.         ErrorHandler(XLM_NOT_SQUARE);
  385.         return (LPFP)NULL;
  386.     }
  387.  
  388.     // init XLK 1 row larger to store eigenvalues in last row
  389.     if ((lpXLKInt =InitRetBuff(wRows+1, wRows)) == NULL)
  390.         return (LPFP)NULL;
  391.  
  392.     // create pointers for XLK internal array which is
  393.     // used to store vectors + eigenvalues in bottom row
  394.     if ((lplpEivec = InitPointers(&lpXLKInt->Data[0], wRows+1, wRows))
  395.             == NULL)
  396.     {
  397.         Success = FALSE;
  398.         goto Error3;
  399.     }
  400.     // create pointers for Hmat
  401.     if ((lplpHmat = InitPointers(&lpHmat->Data[0], wRows, wRows))
  402.             == NULL)
  403.     {
  404.         Success = FALSE;
  405.         goto Error2;
  406.     }
  407.  
  408.     // allocate a Global 1-D array for eigenvalues
  409.     if ((lpEival = (LPREAL)GetMem( sizeof(double)*wRows )) == NULL)
  410.     {
  411.         Success = FALSE;
  412.         goto Error1;
  413.     }
  414.  
  415.     // diagonalize the matrix
  416.     nIter = Jacobi(lplpHmat, lplpEivec, lpEival, wRows);
  417.  
  418.     // copy the eigenvalues
  419.     for (id = 0; id < wRows; id++)
  420.     {
  421.         lplpEivec[wRows][id] = lpEival[id];
  422.     }
  423.  
  424.     // free pointers and arrays
  425.     FreeMem((LPVOID)(LPVOID)lpEival);
  426. Error1:
  427.     MemFreePtr((LPVOID)lplpEivec);
  428. Error2:
  429.     MemFreePtr((LPVOID)lplpHmat);
  430. Error3:
  431.     if (!Success)
  432.     {
  433.         return ((LPFP)NULL);
  434.     }
  435.  
  436.     // return pointer to XL map to EXCEL
  437.     return ( (LPFP) lpXLKInt);
  438. }
  439. /********************************************************************
  440.  PolyCurveFit()
  441.  =============
  442.  This function is a generalized least squares curve fitting function.
  443.  It fits a polynomial with linear coefficients to a dependent -
  444.  independent variable set of data
  445.  *******************************************************************/
  446. LPFP PASCAL FAR _export PolyCurveFit(LPFP lpIndVar, LPFP lpDepVar, unsigned int Order)
  447. {
  448. #define    nXLKCols    3    // placed in a NumObs x nXLKCols array
  449.  
  450.     LPFP lpXLKInt;
  451.     LPREAL     lpYest, lpResid, lpPolycoef, lpCoefsig;
  452.     static double    See, Rsqrval, Rval, dferror;
  453.     WORD wRows = lpIndVar->wRows;
  454.     WORD wCols = lpIndVar->wCols;
  455.     WORD NumObs = wRows > wCols ? wRows : wCols;
  456.     WORD wMaxRows = 2*(Order+1) + 4;
  457.     BOOL Success = TRUE;
  458.     WORD i, wIndex;
  459.  
  460.     if (wMaxRows < wRows)
  461.         wMaxRows = wRows;
  462.  
  463.     // init XLK & zero array data
  464.     if ((lpXLKInt = InitRetBuff(wMaxRows, nXLKCols)) == NULL)
  465.     {
  466.         return (LPFP)NULL;
  467.     }
  468.     for (i=0; i<wMaxRows*nXLKCols; i++)
  469.         lpXLKInt->Data[i] = 0.0;
  470.  
  471.     // memory for estimated Y values
  472.     lpYest = (LPREAL)GetMem( NumObs*sizeof(double));
  473.     if (lpYest == NULL)
  474.     {
  475.         Success = FALSE;
  476.         goto Error4;
  477.     }
  478.     // memory for residuals
  479.     if ((lpResid = (LPREAL)GetMem( NumObs*sizeof(double))) == NULL)
  480.     {
  481.         Success = FALSE;
  482.         goto Error3;
  483.     }
  484.  
  485.     //memory for coefficients of fitted polynomial
  486.     if ((lpPolycoef = (LPREAL)GetMem( (Order+1)*sizeof(double))) == NULL)
  487.     {
  488.         Success = FALSE;
  489.         goto Error2;
  490.     }
  491.  
  492.     // memory for standard errors of coefficient estimates
  493.     if ((lpCoefsig = (LPREAL)GetMem( (Order+1)*sizeof(double))) == NULL)
  494.     {
  495.         Success = FALSE;
  496.         goto Error1;
  497.     }
  498.  
  499.     Success  = _PolyCurveFit((LPREAL)&lpIndVar->Data[0],
  500.         (LPREAL)&lpDepVar->Data[0], (int)NumObs, (int)Order, lpPolycoef,
  501.         lpYest, lpResid, (NPREAL)&See, lpCoefsig, &Rsqrval, &Rval, &dferror);
  502.     if (Success == TRUE)
  503.     {
  504.         // remaining variables stuck in XLKInt directly below polycoef
  505.         wIndex = 2*(Order+1);
  506.         lpXLKInt->Data[(wIndex)*nXLKCols + 2] = See;
  507.         lpXLKInt->Data[(wIndex+1)*nXLKCols + 2] = Rsqrval;
  508.         lpXLKInt->Data[(wIndex+2)*nXLKCols + 2] = Rval;
  509.         lpXLKInt->Data[(wIndex+3)*nXLKCols + 2] = dferror;
  510.  
  511.         // copy data into XLKInt
  512.         for (i=0; i < NumObs;  i++)
  513.         {
  514.             wIndex = i*nXLKCols;
  515.             lpXLKInt->Data[wIndex] = lpYest[i];        // Y estimates
  516.             lpXLKInt->Data[wIndex+1] = lpResid[i];     // residuals
  517.             if ( i <= Order )
  518.             {
  519.                 lpXLKInt->Data[wIndex+2] = lpPolycoef[i];    //
  520.                 wIndex = (i+Order+1)*nXLKCols;
  521.                 lpXLKInt->Data[wIndex+2] = lpCoefsig[i];     //
  522.             }
  523.         }
  524.     }
  525.  
  526.     FreeMem((LPVOID)lpCoefsig);
  527. Error1:
  528.     FreeMem((LPVOID)lpPolycoef);
  529. Error2:
  530.     FreeMem((LPVOID)lpResid);
  531. Error3:
  532.     FreeMem((LPVOID)lpYest);
  533. Error4:
  534.     if (!Success)
  535.     {
  536.         return ((LPFP) NULL);
  537.     }
  538.     return (lpXLKInt);
  539. }
  540. #undef nXLKCols
  541.  
  542. /********************************************************************
  543.  CubicSplines()
  544.  =============
  545.  This function will fit a set of cubic spline polynomial equations
  546.  to a discrete set of data.
  547.  *******************************************************************/
  548. LPFP PASCAL FAR _export CubicSplines(LPFP lpIndVar, LPFP lpDepVar)
  549. {
  550. #define    nXLKCols    4    // placed in a wRows x 4 array
  551.  
  552.     LPFP lpXLKInt;
  553.     WORD wRows = lpIndVar->wRows;
  554.     WORD wCols = lpIndVar->wCols;
  555.     BOOL Success = TRUE;
  556.     WORD i, wIndex;
  557.  
  558.     // init XLK & zero array data
  559.     if ((lpXLKInt = InitRetBuff(wRows-1, nXLKCols)) == NULL)
  560.         return (LPFP)NULL;
  561.  
  562.     Success  = _CubicSplines((LPREAL)&lpIndVar->Data[0],
  563.                             (LPREAL)&lpDepVar->Data[0],
  564.                             (int)wRows-1,
  565.                             (LPREAL)&lpXLKInt->Data[0]);
  566.     if (!Success)
  567.     {
  568.         return ((LPFP)NULL);
  569.     }
  570.     return (lpXLKInt);
  571. }
  572. #undef nXLKCols
  573.  
  574. /********************************************************************
  575.  CalcSpline()
  576.  =============
  577.  This function calculate the cubic spline interpolation of an
  578.  y-value given an x-value and the cubic spline coefficients.
  579.  
  580.  *******************************************************************/
  581. LPREAL PASCAL FAR _export CalcSpline(LPFP lpIndVar, LPFP lpCoef, LPREAL lpX)
  582. {
  583.     LPREAL lpXLEInt;
  584.  
  585.     if ((lpXLEInt = (LPREAL)InitRetBuff(1,1)) == NULL)
  586.         return (NULL);
  587.  
  588.     _CalcSpline((LPREAL)&lpIndVar->Data[0],
  589.                 (LPREAL)&lpCoef->Data[0],
  590.                 (int)(lpIndVar->wRows  - 1),
  591.                 (LPREAL)lpX,
  592.                 (LPREAL)lpXLEInt);
  593.  
  594.     return (lpXLEInt);
  595.  
  596. }
  597. /********************************************************************
  598.  SmoothSG()
  599.  =============
  600.  This function uses the Savitsky - Golay algorithm to reduce the
  601.  noise in a sampled data set.
  602.   *******************************************************************/
  603. LPFP PASCAL FAR _export SmoothSG(LPFP lpData, unsigned int wSmoothNum,
  604.     unsigned int wDerivNum)
  605. {
  606.     LPFP lpXLKInt;
  607.     WORD wRows = lpData->wRows;
  608.     WORD wCols = lpData->wCols;
  609.     WORD wNumDat;
  610.  
  611.     // init XLK & zero array data
  612.     if ((lpXLKInt = InitRetBuff(wRows, wCols)) == NULL)
  613.         return (LPFP)NULL;
  614.  
  615.     // check limits
  616.     if (wSmoothNum < 1 || wSmoothNum > 5)
  617.         return (LPFP)NULL;
  618.     if (wDerivNum > 2)
  619.         return (LPFP)NULL;
  620.  
  621.     // set NumDat
  622.     wNumDat = (wRows > wCols) ? wRows: wCols;
  623.  
  624.     // call smoothing routine
  625.     _DataSmoothSg((LPREAL)&lpData->Data[0],
  626.                     (int)wNumDat,
  627.                     (int)wSmoothNum,
  628.                     (int)wDerivNum,
  629.                     (LPREAL)&lpXLKInt->Data[0]);
  630.  
  631.     return (lpXLKInt);
  632. }
  633. /********************************************************************
  634.  SmoothWeights()
  635.  =============
  636.  This function smooths data via a weighted average
  637.  *******************************************************************/
  638. LPFP PASCAL FAR _export SmoothWeights(LPFP lpData, LPFP lpWeights, LPREAL lpDivisor)
  639. {
  640.     LPFP lpXLKInt;
  641.     WORD wRows = lpData->wRows;
  642.     WORD wCols = lpData->wCols;
  643.     WORD wNumDat;
  644.     WORD wSmoothNum;
  645.  
  646.     // init XLK & zero array data
  647.     if ((lpXLKInt = InitRetBuff(wRows, wCols)) == NULL)
  648.         return (LPFP)NULL;
  649.  
  650.     wNumDat = (wRows > wCols) ? wRows : wCols;
  651.     wSmoothNum = (lpWeights->wRows > lpWeights->wCols) ? lpWeights->wRows :
  652.                     lpWeights->wCols;
  653.  
  654.     _DataSmoothWeights((LPREAL)&lpData->Data[0],
  655.                     (int)wNumDat,
  656.                     (int)wSmoothNum,
  657.                     (LPREAL)&lpWeights->Data[0],
  658.                     (LPREAL) lpDivisor,
  659.                     (LPREAL)&lpXLKInt->Data[0]);
  660.  
  661.     return (lpXLKInt);
  662. }
  663.